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ABSTRACT 

For the first time ever, we derive equations governing the time-evolution of 
fully relativistic slim accretion disks in the Kerr metric, and numerically con- 
struct their detailed non-stationary models. We discuss apphcations of these 
general results to a possible limit-cycle behavior of thermally unstable disks. 
Our equations and numerical method are applicable in a wide class of possible 
viscosity prescriptions, but in this paper we use a diffusive form of the "standard 
alpha prescription" that assumes the viscous torque is proportional to the total 
pressure. In this particular case, we find that the parameters which dominate 
the limit-cycle properties are the mass-supply rate and the value of the alpha- 
viscosity parameter. Although the duration of the cycle (or the outburst) does 
not exhibit any clear dependence on the black hole spin, the maximal outburst 
luminosity (in the Eddington units) is positively correlated with the spin value. 
We suggest a simple method for a rough estimate of the black hole spin based 
on the maximal luminosity and the ratio of outburst to cycle durations. We also 
discuss a temperature-luminosity relation for the Kerr black hole accretion discs 
limit-cycle. Based on these results we discuss the limit-cycle behavior observed in 
microquasar GRS 1915+105. We also extend this study to several non-standard 
viscosity prescriptions, including a "delayed heating" prescription recently stim- 
ulated by the recent MHD simulations of accretion disks. 

Subject headings: accretion, accretion disks — black hole physics — hydrodynamics 
— instabilities 
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1. Introduction 

This work is motivated by tlie fundamental issue of how to find the proper viscosity 
prescription in the black hole (BH) accretion disk theory. The " standard alpha-prescription", 



introduced in the seminal paper by 



Shakura &: SunyaevI (Il973l ). assumes that the viscous 



torque T is proportional to the total pressure P, 

T = -aP. 



where the dimensionless "viscosity coefficient" < a < f is a universal phenomenological 
constant. Hydrodynamical accretion disk models constructed with the help of the above 
formula are remarkably successful in describing observed properties of real disks, but only 
when these disks are stationary and not very luminous, i.e. when 

^ = and L<0.3LEdd. (2) 
This is not a surprise. An extension o f the Shakura-Sunyaev theory in terms of the "slim 



disk" models (lAbramowicz et al 



19881 ). predicted that the properties of accretion disks 
weakly depend on a only in the range ([2]). Outside this range, the predicted observables 
strongly depend on the form of the viscosity prescription and on the assumed value of a: 
the "small-alpha disks", a ~ 0.01, are remarkably different from the "moderate-alpha disks", 
a ~ 0.1. For luminosities higher than 0.3LEdd the dependence on a is convoluted with the 
details of radiative transfer, vertical structure, and general relativistic effects inside and 
outside the disk. In addition, the relativistic effects strongly depend on the BH spin. 

These severe technical complications sometimes mask a fundamental problem. On 
one hand, with the help of hydrodynamical models we are able to accurately calculate 
observable properties of accretion disks, but only if we assume ([T]) or other ad hoc formula 
for the stress. However, outside the range of applicability ([2]), this aproach is questionable 
and most probably wrong. On the other hand, the MHD simulations which calculate T 
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from the first principles, are not sufficiently realistic today in treating radiative transfer 
and other relevant physics accurately enough to calculate realistic observable properties of 
accretion disks. 

One possible way to make progress is to calculate "the state-of-art" hydrodynamical 
models with radiative transfer, vertical structure, general relativistic effects, etc. treated as 
accurately as possible and with a properly chosen "MHD viscosity prescription" motivated 
directly by the MHD simulations. Some phenomenological parameters that may appear in 
the MHD prescription, should be th en fixed by a det ailed comparison wit h obse rvations. 



This principle was recently used by 



Sadowskil (120091 ) and 



Sqdowski et al 



torn who 



calculated a network of models of stationary slim accretion disks, covering a wide range 
of the relevant parameter space. Our present work depends on these models in several 
important aspects. In particular, it takes the stationary models calculated by Sqdowski as 
the initial condition for the non-stationary calculations. 

We also b orrow heavily from the previous paper of our long-term research program 



(ILi et al. 



20071 hereafter Paper I). We have introduced there the adaptive pseudo spectral 



doma in decomposition method, which we also use in the present paper. Results of 



Li et al. 



( 120071 ) have been in a good agreement with th e studies on the limit 



the radiation-pressure supported slim disks by 



SMOl). These, and al 



Szuszkiewicz fc Millerl (12001 



c ycle b ehavior in 



1991 



Lasota and Pelat 



Mayer &: Pringle 



othe r previous numerical studies of th e sub 

] 



1991 



Szuszkiewicz &: Miller 



1997 



1998 



2001 



hereafter 



ect (iHonma et al 



Teresi et al 



2004al Jbl: 



20061 ) . have been performed within the Newtonian hydrodynamics, in 



which the Newtonian gr avitational potentia. 



potential introduced by 



i s repl aced by the Paczyhski pseudo-Newtonian 



Paczyhski fc Wiital (119801 ). It reasonably models gravity of 



non-rotating black holes, but fails in the case of rapidly rotating ones. 



Our present paper extends these previous studies by making three new developments: 
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(i) We perform the first ever fully relativistic time- dependent numerical study for a 
limit-cycle in black hole accretion disks, (ii) We provide a theoretical tool for understanding 
some observational aspects of GRS 1915+105 (and possibly similar sources) which are 
unique in showing both the evidence for a near-extreme black hole spin, and the limit-cycle 
behavior, (iii) We investigate the possibility of using the observed properties of the limit 
cycle to estimate the black hole spin. These issues are relevant for studying the fundamental 
issue of the viscosity prescription. 



In this paper, we consider the axisymmetric relativistic accretion flows around Kerr 
BHs. We use the Boyer-Lindquist spherical coordinates t, r, 6, (j) to describe the space-time 
around a BH. Putting all of the complicated derivations into Appendix \^ the basic 
equations which govern the dynamical behavior of accretion flows can be written as 
following!^ . 

• Mass conservation: 



where S (= 2Hp) is the surface density with H being the half thickness and p being 
the mass density, V is the radial velocity measured in the corotating frame (CRF), 
7 is the Lorentz factor, A = — 2Mr + and A = r'^ + r^c? + IMro? with M 
and a being the BH mass and spin per unit mass, respectively. The derivative of the 



function of ^ and ^ (see the following equations and ^ for the deflnitions of 



2. Basic equations 



'dt 




(3) 



contravariant t-component of the four-velocity ^ is deflned in equation (1A14I) as a 



*Throughout this paper we use the c = G = 1 units. 
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these two derivatives). 

Radial momentum conservation: 

dv_ _ yi - V^A 

'dt ~ 7^1/2 



V dV A 1-V^dp 
+ — 



1 — dr r p dr 



(4) 



where A is defined in equation (lA19ap . p (= -^^pT + Prad) is the total pressure 



consisting of the gas and radiation pressure. 
Angular momentum conservation: 



+ 



dt 7VI - V^A^^ dr ' dr\ dr ) ' ^^"^ 

where C{=u^) is the angular momentum per unit mass, Q is the angular velocity with 
respect to the stationary observer (see equation (lA19bp for its definition), and u is the 



kinematic viscosity coefficient. We evaluate u with a diffusive form of the "standard 
alpha prescription" (equation ([1])), 

u = \aH.I^-. (6) 
3 V P 

Half thickness evolution: 

— = -f/cose,---^=_ — , (7) 

where H is the half thickness of the disk, U is the vertical velocity, and Qh 
(= arccos is the colatitude angle corresponding to the disk surface. 

Surface vertical motion: 

dU ^ U f V dV Cr^dC\ U dH 



dt j^A^/^cosQh -f^\{l-V^ydt A dt J H dt ' 
where TZ is defined in equation (lA33ap . Equations ([7j) and ([8]) describe the 



evolution of the disk surface. This kind of vertical treatment was first introduced 
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by ISzuszkiewicz &: Millerl (l200ll ) in order to mimic so me advantageous features 



of a two-dimensional model in a one-dimensional (see 



Papaloizou &: Szuszkiewicz 



I994J ) study. Thus, our study also exhibits some two-dimensional features (e.g., the 
vertical component of the velocity at the disk surface) so that we prefer to call it 
1.5-dimensional. 



Energy conservation: 



dT 1 



dt S 7AV2 



Cv 



+ (Fs - 1)TS ( 



du^ 1 d 



\ dt dr 



dT 



7VI - dr ' 

(9) 



where T is the temperature of gas, is the local viscous heat generation rate 
(equation (1A37I) ). and F~ is the radiative cooling rate given by the bridge formulae 
which is valid for both optically thick and thin regimes. 



3rR/2 + y3 + l/rp' 

where Tr and Tp are th e Rosseland and Planck mean optical depths (see e.g.. 



(10) 



Abramowicz et al. 



19961). The corresponding radiation pressure is given by. 



Cv and Fa are defined in equations (lA43aP and (lA43bp . respectively. 



:ii^ 



In this paper, the properties of an accretion disk depend on three crucial parameters: 



Dimensionless black hole spin a*: 



a = 



GM/c M' 
Diffusive viscosity parameter a (see equation 



1< a* < 1; 



(12) 



< a < 1; 
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• Dimensionless mass-supply rate rh: 

m = — : . 13 

MEdd 

Where M(rout) and M^dd = Q4ttGM/{ck,cs) are the accretion rate at the outer boundary 
f = ^out and the Eddington accretion rate (corresponding to the Eddington luminosity for a 
non-rotating BH), respectively. 



3. Numerical method 

Following Paper I, we use the adaptive pseudo spectral domain decomposition method 
to solve the basic equations. This method has already been validated by the successful 
reproduction of the limit-cycle behavior in the pseudo-Newtonian framework. The 
relativistic time-dependent code used in this work is an upgraded version of the previous 
pseudo- Newtonian one. For the initia l cond itions, we interpolate the relativistic stationary 
global slim disk solution of Sq.dowski ( 2OO9I ) on the grid of the time-dependent code. We 



study the accretion disk around a BH with mass M = lOM© and setup the computation 
domain from a radius between the BH horizon and ISCO (Innermost Stable Circular Orbit), 
say r = 2.5rg for a nonrotating BH (risco = 3rg), to lO'^rg (rg = 2GM/c^ = 2M). As in 
Paper I, the computational domain includes a highly supersonic region, with high radial 
velocities (about 0.2 ~ 0.3c at the inner boundary). At this inner boundary, we adopt 
the free-type boundary conditions, i.e., we let all variables evolve naturally according to 
their basic equations except C, which is set by an extrapolation of the nearest two points. 
This kind of inner boundary conditions, in practice, are quite stable because the accretion 
flow always supersonically inflows beyond the inner boundary without numerical spurious 
reflection. Meanwhile, they are also more natural and physical than directly setting the 
viscous torque to vanish outside the BH's horizon. At the outer boundary, we fix V and 
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S to ensure constant mass supply, and also fix L and set f/ = to avoid the numerical 
spurious evolution. 

We have performed computations for two particular cases to check the relativistic 
code. Both consist of an accretion disk with a = 0.07 and m = 0.02, while the BH spin 
was set to a* = and a* = 0.95, respectively. The essential difference betwee n these two 



cases is in the thermal stability. According the local analysis of stability (e.g. 



Kato et al 



1998 



p. 155), only the non-spinning case is thermally stable (see Fig. [T]). In Figure |2l we 
show the profiles of surface density and temperature for both cases. It is clear that for 
the a* = case (left panel), the disk remains in the steady state for more than 1.2 x 10^ 
seconds, therefore it can be regarded as a stable disk (the tiny differences between the 
initial state and the snapshot are due to the differences of numerical methods used for 
obtaining stationary and time-dependent solutions and are considered not significant). On 
the contrary, for the highly spinning case (right panel) the disk undergoes an outburst 
immediately after the beginning of simulation. These results validate the ability of the code 
to distinguish the thermally stable and unstable behavior, and also prove that the BH spin 
significantly infiuences the thermal stability of the accretion disk, which motivates us to 
explore the limit-cycle behavior of thermally unstable disks around Kerr BHs. 



4. Numerical exploration 

4.1. Parameter space and characteristic quantities 

In order to understand the effects of BH spin on the limit-cycle behavior, it is necessary 
to investigate the parameter space {a*,a,m) with a series of models. In this section, we 
describe and discuss results of twelve runs assuming different sets of the input parameters 
span on the grid: a* = (0,0.5,0.95), a = (0.07,0.1), and m = (0.06,0.1) (cf. Table [I]). All 
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of the cases are thermally unstable and undergo recursive limit-cycle evolution. 

To describe the limit-cycle behavior, we define four characteristic quantities, which 
are the Outburst Duration (the full width at half maximum of light-curve, hereafter OD), 
Cycle Duration (time interval between two outbursts, hereafter CD), ratio of OD to CD, 
and Maximal Luminosity (the peak value of the light curve during outburst, hereafter ML; 
calculated intrinsically, with no ray-tracing), respectively. In Figure [3l we visualize these 
definitions on a sketched light-curve. 

The computations for each one of those twelve cases continue until the CD converges 
to a constant value. The constancy of CDs implies that the computations has reached a 
state uniquely determined by the parameters (a*,a,m), without any impact of the initial 
conditions. In practice, the required level of convergence is achieved after three or five 
cycles. We have performed computations for each case for more than seven cycles (some 
cases reached ten cycles). Our analysis below is based on the last four cycles only. 

In Table [T] we list the mean values, standard and relative deviations of each 
characteristic quantity for our twelve models. These values are calculated based on the 
retained cycles only. Since the relative deviations (the percentages in the brackets) are all 
very low, we do not show the error bars for the data points in the subsequent figures. 

4.2. Impact of BH spin 

In order to reveal the effects of BH spin on limit-cycle behavior, we have divided our 
twelve models into four groups to refiect the different effects of viscosity and mass supply 
(i.e., mass accretion rate at the outer boundary). In Figures S] and 13 the models with 
(a = 0.07, m = 0.06) and different spins are represented by empty stars, while those with 
(a = 0.1, m = 0.06), (a = 0.07, m = 0.1), and {a = 0.1, m = 0.1) by filled stars, empty 
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circles, and filled circles, respectively (both are based on the data from Tabled]). 

In the four panels of Figure H] we plot OD, CD, OD/CD and ML as a function of a*. 
It is clear that there is no distinct and consistent dependence on a* neither for OD nor CD 
(panels (a) and (b)). The ratio of OD to CD (panel (c)) is almost independent of a* for all 
sets of input parameters. On the contrary, ML exhibits a perfect correlation with a* (panel 
(d)), though this relation is still slightly distorted by the impact of a and m. 

These facts may be easily understood in the framework of the general relativity. The 
effects of BH spin are restricted to the innermost region of accretion disks. Therefore, the 
quantity ML, which corresponds mostly to the radiation emerging from the inner part of an 
accretion disk, must display a strong dependence on a* . The other quantities (OD, CD and 
OD/CD) don't have such a strong dependence on a* as they depend on the disk structure 
and its evolution at a wide range of radii (from the inner part up to more than lOOrg), 
where the effects of viscosity and mass-supply dominates over the effect of BH spin. 

4.3. Possibility of estimating the black hole spin with limit-cycle 

As has been discussed in the previous section, the effects of BH spin on the quantities 
OD and CD are easily disturbed by the complex effects of viscosity and mass-supply; while 
the ratio OD/CD has no discernible dependence on BH spin. Thus, they cannot serve as 
proper probers of BH spin. Even ML, which has monotonic and positive correlation with 
a*, cannot be used directly to estimate the BH spin, because the dependence of ML on BH 
spin is significantly affected by other factors (see panel (d) of Figure Hj). 

None of those quantities can be directly used as a probe, but there still exists the 
possibility of probing the BH spin with them. In fact, we find that the combination of ML 
and OD/CD can be used to estimate the BH spin. In Figured], we show the OD/CD-ML 
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diagram for all the models. The ratio OD/CD raises with increasing values of a and m, and 
is hardly dependent on a* (see panel (c) of Figure H]). It is remarkable that the combined 
effects of viscosity and mass-supply make ML scale almost linearly with the ratio OD/CD. 
For a given BH spin there is a single straight line reflecting the dependence of ML on 
OD/CD. These lines do not intercept and have similar slopes. Therefore, the OD/CD-ML 
diagram may be used for a rough estimation of BH spin if only the limit cycle parameters 
can be obtained from observational data. 

We find that the dependence of ML on OD / CD may be approximated by the following 
formula 

ML / LEdd = 7.59 OD/CD + 0.71 + 7.18 r^, (14) 
where r/ is the efficiency of accretion for a thin disk and is given by 



2M 

r^=l-Wl 



3?"ms 

with Tins being the radius of the marginally stable orbit. In Figure [5] we plot with dashed 
lines the fits obtained with these formulae for a* = 0, 0.5 and 0.95. Given the values of ML 
and OD/CD, one can easily obtain the radius of the marginally stable orbit from Eq. [HI 
and subsequently the BH spin. 

However, there are at least a few limitations for the application of this method. The 
first one is related to the fact that the OD/CD-ML relation is still subject to a significant 
dispersion (see Fig. [5]). For a given BH spin, the combined impacts of viscosity and 
accretion rate result in a relation which is only roughly linear. 

The second one is related to the modulation of the emitted radiation by the 



gravitational redshift and focusing effects near the BH (e.g. ICunningham fc Bardeenlll973l ). 
In our method, the ratio OD/CD is measured in the Boyer-Lindquist (observer's) time, but 
ML is intrinsic, i.e., it is the integral of local radiation fiux in the whole disk. Thus, one 
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should calculate ML only after calculating the observed disk spectrum, which is subject 
to the gravitati onal effects mention ed above. It would be possible by applying ray-tracing 
techniques (e.g. iFanton et al. 



19971 ) to solve this issue. However, the calculations will be 
very complicated and out of the scope of our paper. This is because the precise calculation 
of the emitting spectrum would require knowledge about the intrinsic local emission which 
cannot be approximated as a perfect black body, especially in effectively optically thin 
regions. Even so, we can still anticipate that the apparent MLs would be lower than the 
intrinsic ones (even for the face on case) due to the reduction of the apparent luminosity of 
disks by the gravitational redshift. However, this decrease should not be strong, as most 
of the flux contributing to ML comes from a region extending up to ~ lOr^ where the 
gravitational reddening has little impact. We argue that the observed value of ML is not 
expected to be lower than the intrinsic one by more than a few percent. The observed 
values of ML would be further decreased if the inclination is not face on — one should 
account for this fact before applying the method. 

The last one is the model-dependence of our method. We assume the limited classical 
theory to describe the accretion flow. We neglect mass outflows, magnetic flelds and apply 
the a prescription of viscosity with a independent of radius. These factors, especially the 
viscosity treatment, may have sig niflcant impact on ML vs OD/CD relation that we have 
obtained. iNayakshin et al.l (|2000[ ) show that reproducing the light curves of GRS 1915+105 
is possible only if much more complicated models are considered. Nevertheless, if the 
limit-cycle behavior is the actual reason of the observed variability of some microquasars, 
and if the extraction of ML and OD/CD parameters from the light curves is possible, then, 
under the limitations of our model, the BH spin can be roughly estimated, before other 
more accurate methods, such as fltting spectral energy distribution, are applied. 
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4.4. Evolution of the temperature-luminosity relation for the limit-cycle in 

classical theory 

The maximal disk temperature-luminosity (T-L) correlation predic ted by the classical 



theory of limit-cycles has been us e d to compare with observations (e.g. 



Gierlihski &: Done 



2004 



Kubota fc Makishima 



2004; 



Kubota et al. 



200ll ) . In Figure [6l we show the T-L 



evolution for the last cycle calculated in each model. The panels present the evolution for a 
given set of a and m and different BH spins (a* = 0, 0.5 and 0.95 marked by solid, dashed 
and dotted lines, respectively). These curves reveal several common features: (1) All of the 
cycles begin at the closely-located triangle points and spend a quarter of CD evolving along 
the curve in the anti-clockwise direction before they reach the square points, and finally 
spend the other 3/4 of CD to return to the starting location (triangles); (2) The evolution 
of all the cycles fits (with exception of the outbursts of the a* = 0.95 case) between two 
straight lines correspo nding to Ln^sk oc T^^^^ re lations with different color correction factors 



col 



1 and 7.2, see iGierlinski fc Done 



20041 ): (3) The T-L curves of all the cycles are 



similar for most of CD except for the outburst state — the BH spin impact is revealed at 
the outburst of limit-cycle only; (4) During the mass restoring process (prior to the next 
outburst) the disk obeys Ljisk oc T^'J^ (see the dashed straight lines in all panels). 



5. Discussion 

The theory based on the standard viscosity prescription ([1]) predicts that radiation 
pressure-supported regions of radiatively efficient accretion disks are thermally and viscously 
unstable. The range of luminosity within which the disk is unstable is narrower for an (ad 
hoc) alternative viscosity prescription 

T = -a/P^=^?^, (15) 
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where < /i < 1 is a constant parameter ( jSzuszkiewicz 



19901), with fi = corresponding to 



the standard prescription and /i = 1 to the so-called "geometric al mean" prescriptio n. The 



instability disappears for fi > 8/7. A recent important work by 



Hirose et al. 



(120091 ) shows 



that the standard Shakura-Sunyaev viscosity prescription ([T]) fits their MHD s imulations o: 



accretion disks much better than the prescriptions (fTSl) with /i > 0. However, 



Hirose et al 



( 20091 ) also find that the radiation pressure supported MHD disks are thermally stable in 



their simulations (it is not known whether they are viscously stable). This result contradicts 
the hydrodynamical stabilit y analysis. What ar e the MHD effects that stabilize the disk? 



One explanation is given in 



Hirose et al. 



( 120091 ). They argue that although the viscosity 



prescription has the form ([T]), it should nevertheless be modified in a subtle way, because 
the "viscous heating" that occurs in the disk is delayed with respect to the MHD stress 
T = Tmhd, 



Q+(t) =T(t- At) 



dQ 
dr 



(16) 



Here At is the delay, found to be about 10 dynamical times, and Q is the angular velocity 
of matter. It should be obvious that this modification is irrelevant in the statio nary case. 



In an exhaustiy e Newtonian analytical stability analysis ( ICiesielski et al 



Lin et al. 



(1201 11 ). see also 



(l201l[ )). parallel to the numerical work described in the present pa per, we have 



Hirose et al. 



confir med (in general) that in some parameter space the suggestion made in 
( l2009[ ) does stabilize the radiation pressure thermal instability. However, we also found a 
variety of different parameter ranges with interesting, and complex, oscillatory behaviors 
that need to be further investigated. The present paper makes the first step into this 
direction, by investigating the standard case At = 0. We already numerically simulate a 
few representative cas es with At ^ 0. An other explanation to the absence of the thermal 



instability is given by 



Zheng et al 



(|2011[ ). They attribute the stability to the magnetic 
pressure in the accretion disk, which was neglected in the traditional hydrodynamical 
stability analysis. They show that if the magnetic pressure decreases with response to an 
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increase of temperature of the accretion flow, the threshold of accretion rate above which 
the disk becomes unstable increases significantly compared to the case of not considering 
the magnetic pressure. The physical reason is that in this case the dependence of turbulent 
dissipation heating on temperature becomes weaker. 

On the observational front, spin of BHs in several galactic BH candidates have recently 
been measured by fitting the observed spectral energy distribution with the relat ivistic 



geometrically thin accretion disk model (IShafee et al 



2006 



McClintock et al 



20061 ). Among 



these objects, GRS 1915+105 is the most special one. It is believed to have an e xtreme 



Kerr BH with the spin very close to a* = J/M^ = 1 (see 



McClintock et al 



also t 



1997; 



le only object that exhibits the quasi- regular luminosity variations ( Belloni et al. 



Navaks. 



lin et al. 



Kawata et al 



2006 



2000 



Janiuk et al. 



Janiuk &: Czerny 



2002 



Watarai fc Mineshige 



2003 



2006). It is 



Ohsuga 



2006 



20111), similar to the limit-cycle predicted by the 



standard Newtonian, hydrodynamical models. Several questions need to be investigated 
here. May one explain the luminosity variations observed in GRS 1915+105 as the classic 
"slim-disk" limit-cycle behavior? What is the role of the nearly extreme spin of GRS 
1915+105 in the context of the observed variability? Could the case of GRS 1915+105 help 
to find the answer about the proper viscosity prescription? 

Our research shows that the shape of light curves weakly depend on BH spin. The only 
quantity that feels the impact of the BH spin is the maximal luminosity. Therefore, under 
the assumptions adopted in this paper (e.g., a constant a), the limit cycle can produce light 
curves which are only qualitatively similar to those observed in GRS1915+105, independent 



of BH spin. To o 



Nayakshin et al. 



jtain better agreement one has to construct more complicated models (e.g.. 



20001 ). We prove that the cycle duration weakly depends on the BH spin. 



Thus, the observed value for GRS1915+105 most probably results from combined effects of 
viscosity and mass-supply rate and is not affected by the BH spin. The very presence of 
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li mit cycle variability 



by 



Hirose et al. 



only in this object, if it is true, may challe nge the explanation given 



20091 ) , but could be understood in the model of 



Zheng et al 



mm as due 



to the high luminosity of this source compared to other sources. 



6. Summary and conclusions 



This paper is a sequel of Paper I (ILi et al.l 120071 ). We have changed the framework 



of study from the previous pseudo-Newtonian hydrodynamics to relativistic one. We 
have established the time- dependent basic equations for slim disk accretion in the full 
general relativity (Appendix [A]). We continue to use the adaptive pseudospectral domain 
decomposition method, which has been validated in Paper 1, to solve these new equations. 
Our main results and conclusions are summarized as follows. 

1. We have calculated two models assuming a = 0.07 and m = 0.02 and two values of 
BH spin. We find that the model with a non-spinning BH (a* = 0) can stay in the steady 
state, but the model with a high-spinning BH (a* = 0.95) undergo an outburst immediately 
after the beginning of the computation. These features agree with the predictions of the 
classical local analysis of stability (see Figured]), thus the ability of the new relativistic code 
for distinguishing the thermally stable and unstable disks is validated, and it is confirmed 
that the BH spin can change the thermal stability of an accretion disk. 

2. Following the time-evolution of twelve models we have explored the effects of BH 
spin on the limit-cycle behavior. We define four characteristic quantities to describe the 
limit-cycle behavior: the outburst duration (OD), cycle duration (CD), ratio of OD to 
CD, and maximal luminosity (ML). We find that the effects of BH spin on the quantities 
that depend on the disk structure at a wide range of radii (OD, CD and OD/CD), are 
easily disturbed and even obscured by the combined effects of viscosity and mass-supply. 
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However, ML is an exception. As it depends on the radiation coming from the inner region 
of the disk, it overcomes the impact of viscosity and mass-supply and reveals a monotonous 
positive dependence on BH spin. 

3. We have discussed the possibility of using the OD/CD-ML diagram to estimate 
the BH spin. The advantage of this method is its simplicity, but the dispersion related 
to the unknown viscosity and mass-supply, as well as its model dependence limit its 
application. We suggest that one can use the OD / CD-ML relation to perform rapid and 
rough estimation, before more accurate BH spin estimations based on other methods are 
available. 

4. We have presented and discussed the evolution of the T-L relation for limit-cycles 
with different model parameters, proving that BH spin changes the disk evolution only 
during the outburst phase. 

We thank Feng Yuan for beneficial discussion. This work was supported by 973 
Program under grant 2009CB824800, the National Natural Science Foundation of China 
under grants 10833002 and 11003016, the Natural Science Foundation of Fujian Province of 
China under grant 2010J01017, and by Polish Ministry of Science grants N203 0093/1466, 
N203 304035, N203 380336 and N203 381436. 



A. Derivation of basic equations 



The basic equations used in our paper, which make exploring the ti me-dependent be 



havior of acc r etion disks around K e rr BH s pos sible, are derived basing on 



Abramowicz et al. 



(11996 



19971 ) 



Gammie &: Pophairu ( 119981 ) and 



Sqdowskil (j2009[ ). In this appendix, we give 



all of the derivation details that are ignored in the main part of our paper. 
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A.l. Metric and four-velocity 

We use the Boyer-Lindquist spherical coordinates t, r, 9, (f) to describe the Kerr BH 
metric. We use its reduced form near the equatorial plane, because the disk equations 
involve vertically averaged quantities and the disks we deal with are not very thick. The 
nonvanishing covariant and contravariant components of the reduced metric can be written 
as, 

r-2M 2Ma A ^ , 

9u = — , gtci^ = —, = 9rr = 9ee = r ; (Al) 



and 



,,,_ r-2M A 1_ 



^2 



Where M and a are the mass and specific angular momentum of the BH respectively, 
A = — 2Mr + a^, A = + r^a^ + 2Mra?. The metric is in the geometrical units 
(c = G = 1) and its signature is ( — h ++). 

The stress-energy tensor reads, 

T'^ = pu'vl' + pg^^ + + u^q' + u'q^, (A3) 

where p, u^, p, s^'^ and are the rest mass density, contravariant components of four-velocity, 
total pressure (the sum of the gas and radiation pressure), viscous stress tensor, and 
radiative energy flux, respectively. 

The general form of the four- velocity is, 

(ej,) + V^^^e\,) + y^ej^^ + V^'^e\,)) , (A4) 

where 7 is the Lorentz factor, y^-'^ are velocities as measured in the Local Non Rotating 
Frame (LNRF) and e(j) are LNRF basis vectors. Near the equatorial plane, they take the 
following forms. 
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y{r) 
Y{<I>) 



1 



V 



U COS 9, 
1 Cr 



(A6a) 
(A6b) 
(A6c) 



and 



9 2Ma d 



He) 



r dr ' 
ld_ 

r 89 ' 
r d 



(A7a) 

(A7b) 
(ATc) 
(A7d) 



where we have followed iBardeen et al. 



( I1972I ) to use the relations flA7aP -f lA7d|) . and we 



have introduced the radial velocity measured in the corotating frame (CRF) V , the vertical 
velocity parameter U and the angular momentum per unit mass C{= u^) to describe 
accretion flows. 

According to equations (lA6p . (1A7P and (lA4p . the contravariant components of the 
four-velocity in terms of V , U and C are 



M = 7 



rAV2' 
V AV2 



u 



Cr^ ^1/2 
; — 
r 



= lucos9. 
r 



(A8a) 
(A8b) 

(A8c) 
(ASd) 



The relevant covariant components are 



(A9a) 
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r V 

Ur = ^,in , =1 (A9h) 

= C, (A9c) 

ug = 'jrU COS 9. (A9d) 

At last, according to equation (IA5I) . the Lorentz factor can be written as, 



A. 2. Mass conservation 



The general form of the continuity equation is 

Vi(pMO = 0. 

After vertical integration it takes the following form 



(All) 



u 



dt 



r Or dt ^ 



(A12) 



where S is the surface density. Substituting and with relations (lA8aP and (lA8bp . we 
have 

Pit. rAV2 r ^ F) / v Ai/2\i 

(A13) 



as 

'dt 



dt r or \ a/1 — V"^ r 



ft t 

where the time derivative ^ can be substituted with the equation 



'dt 



^1/2 1 



V 



dV Cr^ dC 

+ 



(A14) 



rAV2^ [(l-\/2)2 ' A dt\' 
which is obtained by operating time derivative on equation (lA8ap . In the right-hand-side 



of the above equation, there still exist two time derivatives ^ and but they can be 



dC 



dt 



dt 



evaluated before the evaluation of ^ with equations ( 1A18I) and ( lA23p . 
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A.3. Radial momentum conservation 

The radial momentum conservation is expressed as the vanishing of the r-component 
of the divergence of the stress-energy tensor 



V,T*'^ = 0. 



(A15) 



Expanding it 

u'V,u' + -g''% = 
p or 

and then neglecting terms proportional to cos^ 6 we get 



1 rr^P 

p or 



This equation can be reformulated as 



dV _ VI - V^A 
'dt ~ 7v4V2 



V dV A 1-V^dp' 
+ 



1 — V"^ dr r p dr 



where A, Q, Q^, Q^, Cl and R are defined as in 



Abramowicz et al. 



fll996[ ) 



A = 

Q = 



MA 

u"^ 2Mar r^A^/^C 
+ 



A 

Ml/2 



7^3/2 



Q = Q 



R = 



r3/2 ± aMi/2 ■ 
2Mar 



A 



A 



.2^1/2 ' 



(A16) 



(A17) 



(A18) 



(Al9a) 

(Al9b) 

(Al9c) 
(Al9d) 
(Al9e) 



We have neglected the radial components of the viscous stress tensor (s*^) as they are 
negligible for vertically integrated models of accretion disks. 
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A. 4. Angular momentum conservation 



The angular momentum conservation can be written as 



v^ine) = 0, 



where ^'^(= S^^^) is the azimuthal Killing vector. After vertical integration we get 



(A20) 



,du 



,du, 



dt 



dr 



(A21) 



where S\ is t 



le vertically integrated (i, 0) component of the viscous stress tensor. We follow 
Lasotal (119941 ) to assume that the only nonvanishing component of is 

^3/2^1/2^3 



dr 



(A22) 



where u is the kinetic viscosity coefficient and Q is the angular velocity with respect to the 
stationary observer [see definition (lA19bp ]. Finally, introducing £ and V with the equations 
I IKM . ( lASbl) and (lA9a we get 



'dt 



dr 



(A23) 



A. 5. Half thickness evolution 

Let us consider the vertical acceleration of the disk surface in the Zero Angular 
Momentum Observer (ZAMO) system of coordinates. The disk surface is defined in 
cylindrical coordinates as 

z = H{r, t) (A24) 

Therefore, the infinitesimal shift in the vertical direction can be expressed in the following 

way 

dz = + ^dt. (A25) 

dr dt 
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Dividing both sides by dt we get 

dz 



dt dr dt 



(A26) 



Introducing U and with the equations (IA6aP and (IA6bp we obtain 



dH ^ 1 V dH 

= -U cosGh j=== — , 

dt 7 V 1 — or 

where 0^ is the colatitude angle corresponding to the disk surface 
H = r cos Qh) at a certain radius. 



(A27) 



Qh and 



A. 6. Surface vertical motion 



We follow lAbramowicz et al.l ( 119971 ) to assume the form of vertically integrated pressure 



as 



P(r,e,t) = Po(r,t) 



cos^ 9 



cos^ Qh 

Let us consider the 6'-component of the divergence of stress-energy tensor 



(A28) 



VuT^ = 0. 



(A29) 



Expanding and vertically integrating we get 



(A30) 



From equation (1A28|) it follows that 



dP 



2P 



cos Qh 



(A31) 



Taking both relations into account we have 



dU 



V dV Cr^ dC 

+ 



dt 72^1/2 cos 72 -1/2)2 dt A dt J H dt' 



U dH 



(A32) 
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where 



n = 



+ {C -a (utut-l))—-^ u — 



S cos Qh 



UtUt = 
r dug 

dr 



7- 



(7t/r cos 9 



(A33a) 
(A33b) 
(A33c) 



A. 7. The energy conservation 



From the general form of the energy conservation equation 



= 0, 



we obtain, in the non-relativistic fluid approximation, 



(A34) 



P 



Similarily like in 



dt dr 



+ V.(r7or**) + V.(mV + mY) = 0. 



(A35) 



Landau fc Lifshita (119591 ) (their Eq. 49.4), after vertical integration, it can 
be reformulated as 



ST 



F-. 



(A36) 



dt dr 

Where 5* is the entropy per unit mass, is the lo c al vis cous heat generation rate and F~ 



Abramowicz et al. 



is the radiative cooling rate (see 



^ r^ \dr ) ^ 

3rR/2 + 1/rp' 
Termodynamical relations give (for p = Prad + Pgas) 

= cyT ( -— - (Ts-1)-^ 
dr \T dr p dr 



(A37) 
(A38) 

(A39) 
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dS ^fldT Idp 



(A40) 



The appropriate sum of the derivatives of p can be substituted using the continuity equation 



tldp ^Idp 
p at par 



dt 



(A41) 



Taking all together we get 



cv 



where 



E u 



'dt 



.dT 
dr 



- (Fa - 1)TE ( 



1 d 



\ dt dr 



r'^u^) 



cy 



4-3^ P 



r, - 1 ET' 



r,-i 



(4-3^)(7ga.-l) 
12(l-^)(7g,,-l)+^' 



P. 



gas 



P 



Ultimately we get 

dT _ 1 r \F+ - F- 
'dt ~ E 7^1/2 [ ^ 



+ (Fa - l)rE 



du^ 
~dt 



l__d_ 



{rW) 



F+ -F- 



(A42) 



VA 



(A43a) 
(A43b) 
(A43c) 

dT 



7V1 - y2^V2 

(A44) 
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Fig. 1. — Prof iles of ratio of gas pressure to total pressure. According to the local analysis of 
stability (e.g. jKato et al.lll998l p. 155), the disk region with (3 < 0.4 is thermally unstable. 
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Fig. 2. — Profiles of surface density and temperature of two special cases. Dash lines cor- 
respond to the initial state, but solid lines are the relevant snapshots at the certain time. 
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Fig. 3.- 



Sketch for explaining the definitions of characteristic quantities. 
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Fig. 4. — Correlation between the characteristic quantities and a*. 
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0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.055 



OD/CD 



Fig. 5. — Plot of maximal luminosity versus ratio of OD to CD. The numbers near each data 
point are the values of a* . The meanings of markers are the same as those in Figure |H The 
dashed lines are the OD/CD-ML fitted lines for different a* cases. 
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Fig. 6. — Evolution of the last cycle for the groups with different a and rh in the maximal 
disk temperature - luminosity (T-L) plane. All cycles begin at the triangle points, and then 
evolve along the curve in the anti-clockwise direction. The square points are used to denote 
the quarter of cycle durations. The solid straight lines in eac h panel represent Ldisk oc T^^x 



Gierlihski fc Done 



20041 ). They are calculated with the 



for a non-spinning BH (see eq. (3) of 
color temperature correction f^oi = 1 for the left lines and fcoi = 7.2 for the right ones. The 
dashed straight lines denote the Ldisk oc T^l^ relation. In each panel, the solid, dashed, and 
dotted curves are for the different BH spins with a* = 0, 0.5, 0.95, respectively. 
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Table 1. Four limit-cycle characteristic quantities of the twelve cases^ 



a* 


a 


rh 


OD [s 




CD [s] 


OD/CD 




ML [LEddl 





0.07 


0.06 


31.21±0.12 


(0.38%)2 


1282±3 (0.23%) 


(2.43±0.02)xl0"^ 


(0.61%) 


1.281±0.001 


(0.07%) 


0.5 


0.07 


0.06 


44.76±0.32 


(0.71%) 


1630±4 (0.25%) 


(2.75±0.03)xl0-2 


(0.96%) 


1.568±0.003 


(0.17%) 


0.95 


0.07 


0.06 


35.72±0.14 


(0.38%) 


1276±10 (0.79%) 


(2.80±0.04)xl0"2 


(1.17%) 


2.219±0.018 


(0.79%) 





0.1 


0.06 


33.78±0.12 


(0.35%) 


976±3 (0.31%) 


(3.46±0.03)xl0-2 


(0.66%) 


1.384±0.003 


(0.18%) 


0.5 


0.1 


0.06 


40.67±0.21 


(0.50%) 


1136±5 (0.44%) 


(3.58±0.04)xl0-2 


(0.94%) 


1.609±0.003 


(0.18%) 


0.95 


0.1 


0.06 


54.68±0.45 


(0.81%) 


1517±10 (0.66%) 


(3.60±0.06)xl0-2 


(1.47%) 


2.415±0.008 


(0.32%) 





0.07 


0.1 


44.80±0.15 


(0.33%) 


1246±4 (0.33%) 


(3.40±0.03)xl0-2 


(0.66%) 


1.385±0.002 


(0.14%) 


0.5 


0.07 


0.1 


46.86±0.12 


(0.24%) 


1233±3 (0.25%) 


(3.80±0.02)xl0-2 


(0.49%) 


1.596±0.003 


(0.17%) 


0.95 


0.07 


0.1 


50.87±0.27 


(0.52%) 


1322±3 (0.23%) 


(3.85±0.03)xl0-2 


(0.75%) 


2.336±0.006 


(0.23%) 





0.1 


0.1 


53.73±0.53 


(0.98%) 


1070±4 (0.38%) 


(5.02±0.07)xl0*2 


(1.36%) 


1.456±0.002 


(0.08%) 


0.5 


0.1 


0.1 


55.47±0.70 


(1.25%) 


1086±9 (0.83%) 


(5.11±0.11)xl0-2 


(2.08%) 


1.692±0.003 


(0.17%) 


0.95 


0.1 


0.1 


53.94±0.47 


(0.87%) 


1084±5 (0.47%) 


(4.98±0.07)xl0-2 


(1.34%) 


2.482±0.004 


(0.16%) 



iSee E]for the definitions of OD, CD, and ML. 

^The values showed in the columns, OD, CD, OD/CD and ML are collected as "mean value "it "standard deviation" 
(relative deviation). 



